Automated borehole geology and petrophysics interpretation using image logs

ABSTRACT

A method of characterizing a borehole traversing an earth formation including (a) obtaining an array of data from a formation characterization tool, wherein the data describes a section of the borehole; (b) computing at least one spatial characteristic describing the relative position of pairs of data; (c) assigning said pairs of data to bins based on the spatial characteristic, wherein the size of the bins are selected based on the tool; (d) transforming the data to petrophysical properties of the borehole; (e) calculating the variance of each bin; (f) developing a model of the variances; (g) determining at least one geostatistical parameter using the model; and (h) upscaling the geostatistical parameters to characterize a region of said earth formation. The method may further include generating a heterogeneity index log using the geostatistical model parameters. The method may be implemented using a computer program product for processing and interpreting borehole data.

The present patent application claims priority to U.S. ProvisionalPatent Application No. 60/468,938, filed May 8, 2003, incorporated byreference herein in its entirety.

FIELD OF INVENTION

The present invention relates to a method to quantitatively describe theformation heterogeneity at the resolution scale of the imagery and, moreparticularly, to a method that applies the techniques of geometricstatistics (geostatistics) to borehole imagery measurements.

BACKGROUND OF THE INVENTION

Borehole imagery is a major component of the wireline business (forexample, Schlumberger's FMI™, Formation MicroScanner, OBMI™ Tools), andan increasing part of the logging while drilling business(Schlumberger's GeoVision™) (as described by Tilke, “Imagery in theExploration for Oil and Gas”, published in Castelli et al., ImageDatabases (2002) Wiley, page 608, incorporated by reference herein inits entirety). While these measurements contain abundant data about thesubsurface, it remains a challenge to automatically extract theimportant geological and petrophysical knowledge contained therein. Theprecise and efficient extraction of this knowledge from the imagesincreases their utility, and therefore will increase the demand forthese measurements.

Schlumberger has a long interest in developing techniques forautomatically interpreting and extracting features from boreholeimagery. Many of the most successful techniques in this area have beenimplemented in Schlumberger's BorTex™, PoroSpect™ and BorDip™applications. BorTex™ semi-automatically segments an image into areas ofsimilar texture that correlate with different rock types. BorTex™ canfurther identify “spots” and “patches” from images for quantification ofvugs and connectivity in carbonates. PoroSpect™ on the other handtransforms an FMI™ conductivity image to a porosity image (see below),then affords the opportunity to analyze the porosity distributions in astatistical sense. Finally, BorDip™ analyzes discontinuities in theimages to automatically identify stratigraphic and structural dips, andfractures.

A limitation of these approaches to automated heterogeneity analysisfrom borehole image interpretation is that they involve treating theborehole image as a two-dimensional image to which image processingtechniques are applied.

While there is a long history of academic and industrial approaches tothe automated interpretation and modeling of borehole imagery (e.g.BorTeX™), there is limited academic or applied work on the applicationof geostatistical techniques for this purpose.

The recognition that borehole images can be mapped to petrophysicalproperties is well understood (see Delhomme et al., “Permeability andPorosity Upscaling in the Near-Wellbore Domain: The Contribution ofBorehole Electrical Images”, (1996) SPE 36822 and Newberry et al.,“Analysis of Carbonate Dual Porosity Systems from Borehole ElectricalImages”, (1996) SPE 35158, incorporated by reference herein in theirentireties).

Perhaps the best established transformation is that from resistivity toporosity space using the classic Archie saturation equation (see Fanchiet al., “Shared Earth Modeling” (2002) Elsevier, page 306, incorporatedby reference herein in its entirety):

$\begin{matrix}{S_{w}^{n} = \frac{{aR}_{mf}}{\Phi^{m}R_{xo}}} & (1)\end{matrix}$where S_(w) is the saturation of the wetting phase (0.00–1.00), n is thesaturation exponent, α depends on the tortuosity (0.35–4.78), R_(mf) isthe resistivity of the mud filtrate, Φ is the porosity (0.00–1.00), m isthe cementation exponent (1.14–2.52), and R_(xo) is the resistivity ofthe flushed zone.

This relationship can be rearranged as follows:

$\begin{matrix}{\Phi = \left( \frac{{aR}_{mf}}{S_{w}^{n}R_{xo}} \right)^{\frac{1}{m}}} & (2)\end{matrix}$

The PoroSpect™ application transforms the resistivity image data intoporosity (Φ) using Equation (2) where the remaining Archie's parametersare either input by the user or derived from other logs. A 1.2-inchvertical window is then applied to these data to generate and analyze aporosity distribution histogram for every depth (0.1 inch verticalspacing) (see Newberry et al., “Analysis of Carbonate Dual PorositySystems from Borehole Electrical Images”, (1996) SPE 35158). Thistechnique has proven very powerful in identifying some aspects ofporosity distribution such as dual-porosity. Note however, thistechnique effectively collapses the image data into a histogram,discarding spatial information (other than depth). It does not considerthe three-dimensional geometry of the FMI™ sensors in the borehole.

An approach to analyze the geostatistical variation of porosity as seenin FMS borehole imagery has been suggested previously (see Delhomme etal., “Permeability and Porosity Upscaling in the Near-Wellbore Domain:The Contribution of Borehole Electrical Images”, (1996) SPE 36822). Itis suggested that the geostatistics on a single pad can be analyzed forfine scale heterogeneity, and intermediate scale heterogeneity can bemodeled across pads.

One object of the present invention is to provide a method to modelformation heterogeneity in terms of geological or petrophysicalproperties.

SUMMARY OF THE INVENTION

The present invention describes a novel approach to quantitativelydescribe the formation heterogeneity at the resolution scale of theimagery by applying the techniques of geometric statistics(geostatistics) to borehole imagery measurements (FMIM, FormationMicroScanner, OBMI™, etc.). It is noted that while the examples providedherein are directed to Schlumberger's FMI™ tool, the method may beapplied to nearly all borehole imagery tools.

The geostatistical parameters, which represent the geological andpetrophysical heterogeneity at the scale of the borehole imagemeasurements, are then used to model the heterogeneity at measurementscales representing larger volumes of investigation (e.g. core plug,porosity logs, resistivity logs). These modeled parameters can then beused to describe the uncertainty in formation properties at any scalegiven measurements taken at a particular scale and are presented assingle channel logs.

An immediate application of this technique is found in carbonates wherelateral heterogeneity calls in to question the utility of conventionallogging tools (e.g. porosity) to determine average formation properties.The technique may also be applied to the interpretation of clastics.

The present invention provides an approach to automated borehole imageinterpretation that treats the borehole as a three-dimensional entity inwhich the tool sensors measure a geological or petrophysical property.This property is thereby modeled in three-dimensions, and the techniquesof modern geometrical statistics (geostatistics) can be applied.

A significant benefit of modeling the three dimensional geostatistics ofthe borehole is that geometric and statistical transformations can thenbe applied to the modeled parameters. An immediate application of thisis the ability to “upscale” the modeled parameters. As will be describedbelow in this document, this allows us to ascribe a heterogeneity index(defined below) to some of the other measurements whose values areaveraged over larger volumes of investigation, e.g. core plugs andlogging tools.

The present invention treats borehole images as three-dimensionalmeasurements of the geological and petrophysical properties of the rock.The technique then uses the established techniques of geostatistics (seeClark, Practical Geostatistics, Elsevier Applied Science Publishers,Ltd. (1984), incorporated by reference herein in its entirety) to modelthe three-dimensional geological and petrophysical heterogeneity of therock in the region of the borehole.

The present invention also discusses the presentation of the modeledparameters as single channel logs summarizing the geology andpetrophysics observed by the imagery.

These geostatistical parameters, which represent the statisticalvariability at the scale of the borehole image measurements, are thenused to model the statistical variability at measurement scalesrepresenting larger volumes of investigation (e.g. core plug, porositylogs, resistivity logs). These modeled parameters can then be used todescribe the accuracy of these larger scale measurements in representingthe mean formation properties and differences expected between them dueto their differing volumes of investigation.

Examples are presented where this technique has been applied to acarbonate well, where lateral heterogeneity calls in to question theaccuracy of conventional logging tools (e.g. porosity).

One embodiment of the present invention is a method of characterizing aborehole traversing an earth formation, comprising: (a) obtaining anarray of data from a formation characterization tool, wherein the datadescribes a section of the borehole; (b) computing at least one spatialcharacteristics describing the relative position of pairs of data; (c)assigning the pairs of data to bins based on the spatial characteristic,wherein the size of the bins are selected based on the resolution of thetool; (d) transforming the data to petrophysical properties of theborehole; (e) determining the variance of each bin; (f) developing amodel of the variances; (g) determining the at least one geostatisticalparameter using the model; and (h) upscaling the geostatisticalparameter to characterize a region of the earth formation. The methodmay further include generating a heterogeneity index log using thegeostatistical model parameters. It is preferred that the array of datadescribes a substantially continuous section of borehole. Whilenon-continuous data may be used, this data may result in “empty” binsand gaps in the geostatical model. The coordinates of the data (relativeto the region of investigation) may be determined, such as based on theborehole geometry, tool geometry, and tool orientation which areobtained as part of the measurements. Accordingly, the spatialcharacteristics describing pairs of data may be based on the distancebetween data pairs, the depth of data pairs in the borehole or withinthe region of investigation, and the orientation of the data pairs (i.e.the azimuth). The data may be upscaled to a three-dimensional modelusing the relative or absolute distance between data pairs and theorientation or azimuth of the data pairs.

The variance model may be developed by (a) computing the variance of thespatial characteristic of each bin; (b) computing an experimentalsemi-variogram using the variances; (c) deriving a model semi-variogramfrom the experimental semi-variogram; and (d) determining thegeostatical parameters using the model semi-variogram.

The method of the present invention may be implemented using a computerprogram product for processing and interpreting borehole data,comprising a computer useable medium having computer readable programcode embodied in the medium for processing borehole data, the computerprogram product having (i) computer readable program code means forcomputing spatial characteristics describing the relative position ofpairs of data, wherein the data describes a borehole; (ii) computerprogram code means for assigning said pairs of data to bins based on thespatial characteristic, wherein the size of the bins are selected basedon the tool used to collect the data; (iii) computer program code meansfor transforming the data to petrophysical properties of the borehole;(iv) computer program code means for determining the variance of eachbin; (v) computer program code means for developing a model of thevariances; (vi) computer program code means for determining at least onegeostatistical parameter using the model; and (vii) computer programcode means for upscaling the geostatistical parameter to characterize aregion of the earth formation. The computer program may further comprisea computer program code means for generating a heterogeneity index logusing the geostatistical model parameter.

Further features and applications of the present invention will becomemore readily apparent from the figures and detailed description thatfollows.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1( a) is a flow chart showing one embodiment of the method of thepresent invention; FIG. 1( b) shows one implementation of the presentinvention.

FIG. 2( a) is a schematic of the FMI™ Tool having the four pairs of padsand flaps generating 192 measurements per depth; FIG. 2( b) is a detailof one pad.

FIG. 3 is a plan view of borehole illustrating (a) the distribution ofFMI™ sensors around perimeter of a 8.5 inch diameter borehole and (b)the vector of length L separating sensors S₁ and S₂ with angularseparation θ in hole of radius R is given by

$L = {2\; R\mspace{11mu}{{\sin\left( \frac{\theta}{2} \right)}.}}$

FIG. 4 is a schematic illustration of a three dimensional extension ofthe present invention.

FIG. 5 is a histogram illustrating frequency of sensor pairs per 0.2inch sampling bin for 8.5 inch diameter hole; Bin 1 is 0–0.2 inch, Bin 2is 0.2–0.4 inch, etc.

FIG. 6 is a horizontal semi-variogram (for the middle of the interval ofFIG. 12) illustrating experimental data from the FMI™ as well as themodeled semi-variogram for FMI™, core plug, a high resolution tool, anda low resolution tool.

FIG. 7 is a conceptual diagram of a series of 11 core plugs and theirmeasured porosity values taken from a homogeneous carbonate.

FIG. 8 is a conceptual diagram of a series of 11 core plugs and theirmeasured porosity values taken from a heterogeneous carbonate.

FIG. 9 is a conceptual diagram illustrating the relative volumes of FMI™measurements, core plugs, and a density log.

FIGS. 10( a) and (b) are graphs of modeled semi-variogram parameters asa function of depth and (a) lag (L) and (b) sill (C) for various volumesof investigation ν: FMI™, Core Plug, High-resolution density porosity,Low-resolution density porosity.

FIG. 11 is a log of heterogeneity index Ψ_(Φ,ν) as a function of depthin the borehole.

FIG. 12 is a porosity image of 30 inch interval with fossiliferous layercontaining 1–2 inch diameter rudists and having a mean porosity ofapproximately 18%.

FIG. 13 is a porosity image of 22 inch high interval exhibiting meanporosity of approximately 22%.

FIG. 14 is a horizontal semi-variogram for the middle of the intervaldisplayed in FIG. 13.

FIG. 15 is a composite porosity log displaying logging tool porositywith core plug porosity.

DETAILED DESCRIPTION OF THE INVENTION

As noted above, one embodiment of the present invention provides anoutput of a series of single channel logs describing the geological andpetrophysical heterogeneity of a borehole. While porosity (Φ) will beused as a measure of this heterogeneity, and FMI™ imagery will be usedas the source measurements, one skilled in the art will recognize thatother parameters and imagery sources may be suitably employed inaccordance with the present invention.

In one embodiment, the generation of the heterogeneity logs involves thefollowing steps, as shown in FIG. 1( a):

-   -   1. Obtain an array of data describing (or imaging) a section of        the borehole wall, for example, a 2D array of pixels describing        the borehole wall 10.    -   2. Compute the spatial characteristics describing the relative        position of the pairs of data 15. This may include calculating        the coordinates (such as the 3D Cartesian coordinates) of the        data from hole geometry, tool geometry, and tool orientation to        determine the depth of, distance between or orientation        (azimuth) of the data pairs.    -   3. Identify appropriate bins, such as based on the spatial        characteristics selected and/or the resolution of the tool 20.        This may be performed for a given depth(s). The size of these        bins will be dependent upon the tool used to acquire the data.    -   4. Assign data pairs to bins representing similar spatial        relationships 25.    -   5. Transform the data values as desired 30, for instance        normalized FMI™ resistivity (R_(xo)) pixels may be transformed        to porosity (Φ) pixels. Other possibilities include but are not        limited to grouping the pixels by their measurement value into        bins using either new or preexisting image analysis techniques.    -   6. Determine the Φ variance for each bin 35 and develop a model        of the variance 40. This variance model may be developed by        computing the variance of each bin 40 a, computing an        experimental semi-variogram 40 b by computing the Φ variance for        each bin, then deriving a model semi-variogram from the        experimental semi-variogram 40 c, thereby determining the        modeled geostatistical parameters 40 d and making them available        for further processing or display.    -   7. Determine the geostatistical parameters from the modeled        variance 45 and compute the same parameters for measurements        with larger volumes of investigation using upscaling approaches        50.    -   8. The up-scaled geostatistical model parameters can then be        used to generate heterogeneity index logs for corresponding core        and logging measurements 55.

These methods may be implemented on a computer readable medium, such asone shown in FIG. 1( b) having a program carrier 60, a storage device65, a processing unit 70, a keyboard (75) or other data input mechanismand a display 80 or other output mechanism.

These steps and potential extensions are described in greater detail inthe following sections. While the example provided herein is a verticalwell of constant width having horizontal layers, the method may beeasily extended to deviated wells with varied width, havingnonhorizontal dip and azimuth.

Compute Three-Dimensional Sensor Geometry

Borehole imaging tools produce images of the measured parameter(resistivity, acoustic impedance, density, etc.) on the rock face. At agiven depth, the exact geometry of the pixels is therefore a function ofthe hole diameter at that location and the sensor configuration. Thus,given the orientation of the borehole, the azimuthal orientation of thetool, and the diameter of the hole (all of which are recorded), the 3DCartesian coordinate of every pixel is known.

Note that while the technique presented here is general, for simplicitythe FMI™ Tool is used to illustrate the technique. It is noted, however,that the present methodology may be easily adapted for use with otherimaging tools.

The FMI™ Tool (shown in FIGS. 2( a) and (b)) is comprised of four pairsof pads and flaps. Each pad and flap (see FIG. 2( b)) has 24 resistivitysensors resulting in 192 sensors per depth (8×24). The four arms of thetool expand in the borehole so that the pads and flaps press against therock face. At a given depth, the sensor geometry is therefore a functionof the hole diameter at that location (shown in FIGS. 3( a) and (b)).Thus, given the orientation of the borehole, the azimuthal orientationof the FMI™ Tool, and the diameter of the hole (all of which arerecorded); the 3D position of every sensor is known.

While the examples presented here in are from a vertical well ofconstant diameter (8.5 inches) to simplify the analysis, the methodologymay be extended to other configurations. FIG. 4 is a schematic diagramillustrating a section of a borehole (depicted as a cylinder)intersected by a geological bed (Plane B). The vertical Plane A iscoincident with the borehole axis and the maximum apparent dip (δ) ofthe bedding relative to the borehole. Knowledge of the sensor positionsis useful to extend the theory illustrated in FIG. 3( b) and Equation(5) below on the ellipse formed by the intersection of Plane B with theborehole (shaded area in FIG. 4). The two parameters r (apparent radius)and β (apparent azimuth) are computed. The following equations expressthese values in terms of a (the azimuth of the sensor), δ (dip to thesensor), and R (the radius of the borehole) which may be determinedusing tool measurements:r ² =R ²(1+cos²α tan²δ)  (3)

$\begin{matrix}{{\sin^{2}\;\beta} = \frac{\sin^{2}\alpha}{1 + {\cos^{2}\;\alpha\mspace{11mu}\tan^{2}\;\delta}}} & (4)\end{matrix}$Identify Two-Point Sensor Pairs

For the two-point geostatistical analysis, sensor pairs should beidentified. Measurements are obtained over a continuous logging run;however, for simplicity, this example considers data collected fromsensors at over a determined depth, and ignores the orientation of thevectors between the sensors (only their magnitude is considered). Thefull three-dimensional analysis is a straightforward extension of thisreduced approach.

As noted above, there are 192 sensors per depth in the FMI™ Tool. Assuch, there are

$\frac{192 \times 191}{2} = 18336$possible sensor pairs. Once the 3D coordinate of each sensor is known,then the Euclidean distances between sensor pairs are triviallycalculated (see FIGS. 3( a) and (b)):

$\begin{matrix}{L = {2\; R\mspace{11mu}{\sin\left( \frac{\theta}{2} \right)}}} & (5)\end{matrix}$Assign FMI™ Sensor Pairs to Bins with Similar Lags

For geostatistical analysis, these sensor pairs should be grouped bypairs into bins with similar L (lags, also referred to as relativedistance). The maximum number of bins possible to achieve maximumresolution of the semi-variogram structure is preferred, but too small abin size will result in aliasing (see Bloomfield, “Fourier Analysis ofTime Series: An Introduction”, 2^(nd) Edition (2000) Wiley, page 261,incorporated by reference herein in its entirety). The correct minimumbin size is dictated by the Nyquist Frequency which is twice the sensorspacing (0.1 inches), i.e. 0.2 inches. For the 8.5-inch diameter holeexample, this results in 43 bins of equal size (except the last which is0.1 inches). Because of the trigonometric effects and non-uniformdistribution of sensors as illustrated in FIGS. 3( a) and (b), thenumber of sensors in each of the 43 bins is not uniform as shown in FIG.5.

Transform Normalized FMI™ Resistivity Image to Porosity Image

Once the sensor pairs and bins have been identified, the resistivityimage is mapped to a geological or petrophysical property, so that thisproperty can be analyzed in a geostatistical sense. While it is possibleto use any attribute extracted from the resistivity image, for thepurposes of this document, PoroSpect™ will be mirrored and Archie'stransformation will be applied to analyze the porosity Φ of a each pixelof the image (see Equation (2)).

In applying Archie's transformation to the entire well, it is assumedthat the Archie parameters in Equation (2): S_(w), m, n, α, and R_(mf)are constant over a given depth window in the borehole

This assumption allows Equation (2) to be expressed as:Φ=λR _(xo) ^(−γ) ^(m)   (6)where

$\begin{matrix}{\lambda = \left( \frac{{aR}_{mf}}{S_{w}^{n}} \right)^{\frac{1}{m}}} & (7)\end{matrix}$

One of the assumptions, therefore, is that λ is constant for a givendepth, but it can vary along the length of the well.

Rather than estimating λ from unknown Archie parameters, it can beeliminated from Equation (6) in the following manner.

From Equation (6), and the previous assumptions, at a given depth:

$\begin{matrix}\begin{matrix}{\left\langle \Phi \right\rangle = {\lambda\left\langle R_{xo}^{- \frac{1}{m}} \right\rangle}} \\{{where}\mspace{14mu}\left\langle {\;\;} \right\rangle}\end{matrix} & (8)\end{matrix}$denotes the expectation or mean of the enclosed expression.

Equations (7) and (8) can be combined to eliminate λ yielding:

$\begin{matrix}\begin{matrix}{\frac{\Phi}{\left\langle \Phi \right\rangle} = \frac{R_{xo}^{- \frac{1}{m}}}{\left\langle R_{xo}^{- \frac{1}{m}} \right\rangle}} \\{{where}\mspace{14mu}\left\langle R_{xo}^{- \frac{1}{m}} \right\rangle}\end{matrix} & (9)\end{matrix}$represents the mean value of R_(xo) ^(−γ) ^(m) at the depth of interestand can be obtained from the resistivities of the image pixels, and

⟨Φ⟩represents the mean porosity at the given depth and can be obtained fromconventional low resolution porosity logging tools. This approach ofusing a low resolution measurement to calibrate a high resolutionmeasurement has been described previously (see Flaum et al., “EnhancedResolution Processing of Compensated Neutron Logs”, (1986) SPE 15541(incorporated by reference herein in its entirety).

Note, due to heterogeneity

⟨Φ⟩has its own heterogeneity index (defined below) because it is a loggingmeasurement. This uncertainty can be integrated into the finalheterogeneity calculations.Compute Experimental Semi-Variogram

Now that the porosity data has been acquired and the bins defined,sample variance for each bin h(i) can be computed with the conventionalexperimental variogram equation (see Clark, Practical Geostatistics,Elsevier Applied Science Publishers, Ltd. (1984)):

$\begin{matrix}{{{2\;{\gamma^{*}\left\lbrack {h(i)} \right\rbrack}} = \left\langle \left( {\Phi_{t + h} - \Phi_{t}} \right)^{2} \right\rangle};{i = {\left\{ {1,2,\;\ldots\mspace{11mu},k} \right\}\mspace{14mu}{where}\mspace{14mu}\left\langle \; \right\rangle}}} & (10)\end{matrix}$again denotes expectation or mean, and k is the number of bins.Model Semi-Variogram

For the purposes of this analysis, a simple semi-variogram model hasbeen defined to fit the experimental data from Equation (10):γ(h)=C ⁰ +C[1−e ^(−(h/L)) ² ]  (11)Thus, the model semi-variogram γ is a function C⁰ (nugget effect) andthe Gaussian model C (sill) and L (range or correlation length).

Over the past 20 years that has been a great deal of research into theoptimal way to fit a model semi-variogram of Equation (11) to theexperimental semi-variogram of Equation (10) (see Cressie, “FittingVariogram Models by Weighted Least Squares”, Mathematical Geology (1985)17, pages 563–588 (incorporated by reference herein in its entirety) andZhang et al., “On the Weighted Least-Squares Method for Fitting aSemivariogram Model”, Computers and Geosciences (1995) 21, pages605–608). It was found that the method of Zhang et al. yields the mostsatisfactory results. This method involves minimizing the cost criterionJ (λ):

$\begin{matrix}{{J(\lambda)} = {\sum\limits_{i = 1}^{k}{\frac{N_{h{(i)}}}{{h(i)}^{2}}\left\lbrack {{\gamma^{*}\left( {h(i)} \right)} - {\gamma\left( {h(i)} \right)}} \right\rbrack}^{2}}} & (12)\end{matrix}$where N_(h(i)) is the number of pairs in bin h(i).

In FIG. 6, the resulting model semi-variogram is plotted along with theexperimental semi-variogram. Superimposed on the experimental dataobtained from the porosity image are model variograms for FMI, CorePlug, High Resolution Tool, and Low Resolution Tool. In the experimentaldata, note the sill with Φ²=0.00077 (Φ=2.8%) reached at a lag ofapproximately 2 inches. This analysis can be repeated at every depthwith imagery in the borehole. The resulting vectors of C⁰, C, and L canthen be plotted as conventional logs. Examples of this is illustrated inFIGS. 10( a) and (b).

Upscaling

Now that the experimental semi-variogram observed at the small (e.g.FMI™) scale has been modeled, the model may be upscaled to other volumesof investigation.

FIG. 7 is a conceptual design of a series of 11 core plugs and theirmeasured porosity values taken from a homogeneous carbonate. For thehomogeneous core plugs, the measured porosity at the core plug scaledoes not vary. By contrast, for the heterogeneous core plugs of FIG. 8,the measured porosity at the core plug scale is highly variable,depending on whether an individual plug intersects a large hole or a lowporosity cemented region. A porosity measurement with a larger volume ofinvestigation (such as a density tool) would read an average of a fewplugs, thus having a smaller variance.

In addition to FMI™ and core plug measurement volumes, two logging toolmeasurement volumes have been considered in this analysis: the highresolution and low resolution porosity tools. The modeled dimensions ofthese various measurement volumes are listed in Table 1.

TABLE 1 Dimensions of modeled rectangular prism volumes of variousmeasurements Width Height Depth Volume Volume Measurement (in) (in) (in)(in³) (relative) FMI ™ Sensor 0.1 0.1 0.1 0.001     1x Core Plug 1.0 1.01.0 1.0  1,000x High Resolution Tool 3.0 3.0 3.0 27.0  27,000x LowResolution Tool 3.0 12.0  4.0 144.0 144,000x

FIG. 9 is a conceptual diagram of the relative volumes (scale effects)of FMI™ measurements, core plugs and a density tool. A set of porositymeasurement taken with the density tool at the location of each coreplug reads the average of numerous core plugs and will thereby have asmaller variance than the same set of core plugs. Because the averagemeasurements taken by the porosity tool will cover overlapping regions,there will not be a large variance between measurements. By contrast,FMI™ measurements have a larger variance than core plug measurements.

If two finite volumes of investigation ν (small volume) and V (largevolume) are considered, then it is possible to derive the geostatisticalparameters for V (i.e. C_(V) ⁰, L_(V) and C_(V)) from the modeledgeostatistical parameters for ν (i.e. C_(ν) ⁰, L_(ν) and C_(ν)).

To upscale the nugget effect C⁰ in Equation (11) it can be shown(Frykman et al., “Geostatistical Scaling Laws Applied to Core and LogData”, (1999) SPE 56822, incorporated by reference herein in itsentirety):

$\begin{matrix}{C_{V}^{0} = {C_{v}^{0}\frac{v}{V}}} & (13)\end{matrix}$

It can also be shown (see Frykman et al.) that the range (L) in Equation(11) increases as a function of the increase in volume size:L _(V) =L _(ν)+(|V|−|ν|)  (14)where |V| and |ν| are the dimensions (lengths) of the volumes V and ν inthe direction of L.

Finally, to upscale the modeled sill C in Equation (11), the concept ofthe point scale sill (C_(p)) is introduced. C_(p) is defined as the sillfor an infinitesimally small volume of investigation. Then, for thesmall finite volume ν, the decrease in sill is defined as follows:C _(p) −C _(ν) =C _(p){overscore (Γ)}ν  (15)where C_(p){overscore (Γ)}ν is the point-scale sill within the volume ν,and {overscore (Γ)}ν is the normalized point-scal sill in ν (definedbelow in Equation (17)).

Equations (14) and (15) are then combined to eliminate the point scalesill C_(p) yielding the following definition for C_(V):

$\begin{matrix}{C_{V} = {C_{v}\frac{1 - {\overset{\_}{\Gamma}}_{V}}{1 - {\overset{\_}{\Gamma}}_{v}}}} & (16)\end{matrix}$

The normalized point-scale sill in ν is obtained from the double volumeintegral:

$\begin{matrix}{{\overset{\_}{\Gamma}}_{v} = {\frac{1}{v^{2}}{\int_{v}^{\;}{\int_{v}^{\;}{{\Gamma_{p}\left( {{\overset{\_}{x} - \overset{\_}{y}}} \right)}\ {\mathbb{d}x}\ {\mathbb{d}y}}}}}} & (17)\end{matrix}$where ┌x−y┐, is the Euclidean distance between points {overscore (x)}and {overscore (y)}, and Γ_(p) is the normalized point scale modelsemi-variogram expressed as:

$\begin{matrix}{{\Gamma_{p}\left( {{\overset{\_}{x} - \overset{\_}{y}}} \right)} = {1 - {\mathbb{e}}^{- {(\frac{{\overset{\_}{x} - \overset{\_}{y}}}{L_{v} - {v}})}^{2}}}} & (18)\end{matrix}$where L_(ν)−|ν| is the modeled point scale range form Equation (11).

Similarly, {overscore (Γ)}ν is expressed as:

$\begin{matrix}{{\overset{\_}{\Gamma}}_{V} = {\frac{1}{V^{2}}{\int_{V}^{\;}{\int_{V}^{\;}{{\Gamma_{p}\left( {{\overset{\_}{x} - \overset{\_}{y}}} \right)}\ {\mathbb{d}x}\ {\mathbb{d}y}}}}}} & (19)\end{matrix}$

Having upscaled the geostatistical parameters to the larger measurementvolumes (Table 1) over the entire relevant depth range of the borehole,it is now possible to present these data as well logs (see FIGS. 10( a)and (b)). Note that L increases as ν increases and C decreases as νincreases.

Construct Heterogeneity Index Logs for Porosity Measurements

For the purposes of this discussion, heterogeneity Index Ψ_(V) isdefined as:

The standard deviation due to heterogeneity of a measurement thatsamples a volume V.

In other words, measurements sampling a volume V over volumes separatedby a distance greater than L_(V) (Equation (14)) will vary with astandard deviation Ψ_(V). Further, the sample variance due toheterogeneity is then given by Ψ_(V) ².

For a simple model semi-variogram, the magnitude of the sill isequivalent to the sample variance (see Clark, Practical Geostatistics,Elsevier Applied Science Publishers, Ltd. (1984)). Ψ_(V) can thereforebe expressed in terms of the proposed upscale parameters C_(V) ⁰ andC_(V).Ψ_(V) ² =C _(V) ⁰ +C _(V)  (20)

Thus, Equation (20) allows the computation of the heterogeneity indexfor all of the idealized volumes of Table 1 at each imaged depth in thewell. This is illustrated by the example heterogeneity index log(Ψ_(Φ,V) as function of depth) in FIG. 11. Note that this intervalexhibits significant heterogeneity at all scales; further, Ψ is as highas 12% for the FMI scale and 9% for the core plug scale.

Porosity logs for a particular tool (e.g. density porosity), ormeasurement (e.g. core plug) can then be expressed with heterogeneityindex envelopes as shown in FIG. 15.

DISCUSSION

FIG. 13 is a section of the reservoir having moderate porosity (22%) andrelatively small objects less than an inch in diameter. Thesemi-variogram from this section (see FIG. 14) illustrates a relativelyshort range (approximately 0.4 inches) corresponding to the visual sizeof the objects. [This can be compared to the interval presented above inFIG. 12, with a relatively long range of 0.86 inches.]

While the different variogram ranges appear to be consistent with thedifferent geology in these intervals, the effect on the upscaledheterogeneity is perhaps more interesting. The short range seen in FIG.14 is responsible for the relatively small upscaled heterogeneity indexcompared to that in FIG. 6. Note that the FMI™ scale sill with(Φ²=0.00029(Φ=1.7%) is reached at a lag of approximately 1 inch. Alsonote that the very low upscaled sills, e.g. Φ²_(coreplug)=0.000036(Φ=0.6%), result from the effect of the relativelylow range (as compared to FIG. 6).

The effect is computed from the upscaling equations. This predictedrelationship between short range yielding greater suppression inupscaled heterogeneity also intuitively makes sense: small objects willnot be resolved when the sample volume is large, while large objectswill be resolved. Thus, it is clear that an understanding of thegeostatistics and varying ranges is preferred to reconcile theheterogeneity observed at the different scales of measurement.

A comparison of the predicted heterogeneity at the core plug scale ofmeasurement with actual core plug porosity measurements is illustratedin FIG. 15. In both the heterogeneous and homogeneous intervals, thepredicted and observed core plug heterogeneity are remarkablyconsistent. Therefore, the apparent discrepancy can be explained byheterogeneity. Also shown in this figure is the ±2Ψ_(core plug) envelopearound the logging tool porosity. It is noted that, as predicted,approximately 2σ of the core plug measurements fall within theheterogeneity index envelope. Further, the log straddles a heterogeneouszone above (0–32 feet) and a homogeneous zone below (32–70 feet). Thedifferences in heterogeneity between the two zones are reflected in boththe heterogeneity index log and the core plug measurements.

While the invention has been described herein with reference to certainexamples and embodiments, it will be evident that various modificationsand changes may be made to the embodiments described above withoutdeparting from the scope and spirit of the invention as set forth in theclaims.

1. A method of characterizing a borehole traversing an earth formation,comprising: a) obtaining an array of data from a formationcharacterization tool, wherein said data describes a section of saidborehole; b) computing at least one spatial characteristic describingthe relative position of pairs of data; c) assigning said pairs of datato bins based on said spatial characteristic, wherein the size of saidbins are selected based on said tool; d) transforming said data topetrophysical properties of said borehole; e) calculating the varianceof each bin; f) developing a model of said variances; g) determining atleast one geostatistical parameter using said model; and h) upscalingsaid geostatistical parameters to characterize a region of said earthformation.
 2. The method of claim 1, further comprising: i) generating aheterogeneity index log using said geostatistical model parameters. 3.The method of claim 1, wherein said array of data describes asubstantially continuous section of said borehole.
 4. The method ofclaim 1, wherein said array of data is an array of pixels.
 5. The methodof claim 1, wherein said array of data is selected from the groupconsisting of resistivity data, density data, and acoustic impedancedata.
 6. The method of claim 1, further comprising determining thecoordinates of said data.
 7. The method of claim 6, wherein saidcoordinates are based on the borehole geometry, tool geometry, and toolorientation.
 8. The method of claim 7, wherein said spatialcharacteristics describes at least one of the distance between datapairs, the depth of said data pairs, and the orientation of said datapairs.
 9. The method of claim 1, wherein said spatial characteristicinclude the distance between data pairs and the orientation of said datapairs, and wherein said upscaling further includes developing a threedimensional characterization of said earth formation.
 10. The method ofclaim 1, wherein developing a model of the variances includes: a)computing the variance of the spatial characteristic of each bin; b)computing an experimental semi-variogram using said variances; c)deriving a model semi-variogram from said experimental semi-variogram;and d) determining the geostatical parameters using said modelsemi-variogram.
 11. A computer program product for processing andinterpreting borehole data, comprising: a) a computer useable mediumhaving computer readable program code embodied in said medium forprocessing borehole data, said computer program product having: b)computer readable program code means for computing at least one spatialcharacteristic describing the relative position of pairs of data,wherein said data describes a borehole; c) computer readable programcode means for assigning said pairs of data to bins based on saidspatial characteristics, wherein the size of said bins are selectedbased on the tool used to collect the data; d) computer readable programcode means for transforming said data to petrophysical properties ofsaid borehole; e) computer readable program code means for determiningthe variance of each bin; f) computer readable program code means fordeveloping a model of said variances; g) computer readable program codemeans for determining at least one geostatistical parameter using saidmodel; and h) computer readable program code means for upscaling saidgeostatistical parameters to characterize a region of said earthformation.
 12. The computer program product of claim 11, furthercomprising: i) computer program code means for generating aheterogeneity index log using said geostatistical model parameter. 13.The computer program product of claim 11, wherein said array of datadescribes a substantially continuous section of said borehole.
 14. Thecomputer program product of claim 11, wherein said array of data is anarray of pixels.
 15. The computer program product of claim 11, whereinsaid array of data is selected from the group consisting of resistivitydata, density data, and acoustic impedance data.
 16. The computerprogram product of claim 11, wherein said code means records thecoordinates of said data.
 17. The computer program product of claim 16,wherein said coordinates are based on the borehole geometry, toolgeometry, and tool orientation.
 18. The computer program product ofclaim 17, wherein said spatial characteristics describes at least one ofthe distance between data pairs, the depth of said data pairs, and theorientation of said data pairs.
 19. The computer program product ofclaim 11 wherein said computer readable program code means for upscalingfurther includes means for developing a three dimensionalcharacterization of said earth formation.
 20. The computer programproduct of claim 11, computer readable program code means for developinga model of said variances includes program codes means to: a) computethe variance of the spatial characteristic of each bin; b) compute anexperimental semi-variogram using said variances; c) derive a modelsemi-variogram from said experimental semi-variogram; and d) determinethe geostatical parameters using said model semi-variogram.